This repository has been archived by the owner on Oct 14, 2023. It is now read-only.
/
Propagation using Cowells formulation.mystnb
210 lines (154 loc) · 5.49 KB
/
Propagation using Cowells formulation.mystnb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
---
jupytext:
text_representation:
extension: .mystnb
format_name: myst
format_version: 0.13
jupytext_version: 1.14.0
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
# Cowell's formulation
For cases where we only study the gravitational forces, solving the Kepler's equation is enough to propagate the orbit forward in time. However, when we want to take perturbations that deviate from Keplerian forces into account, we need a more complex method to solve our initial value problem: one of them is **Cowell's formulation**.
In this formulation we write the two body differential equation separating the Keplerian and the perturbation accelerations:
$$\ddot{\mathbb{r}} = -\frac{\mu}{|\mathbb{r}|^3} \mathbb{r} + \mathbb{a}_d$$
+++
<div class="alert alert-info">For an in-depth exploration of this topic, still to be integrated in poliastro, check out [this Master thesis](https://github.com/Juanlu001/pfc-uc3m).</div>
+++
<div class="alert alert-info">An earlier version of this notebook allowed for more flexibility and interactivity, but was considerably more complex. Future versions of poliastro and plotly might bring back part of that functionality, depending on user feedback. You can still download the older version [here](https://github.com/poliastro/poliastro/blob/0.8.x/docs/source/examples/Propagation%20using%20Cowell's%20formulation.ipynb).</div>
+++
## First example
Let's setup a very simple example with constant acceleration to visualize the effects on the orbit:
```{code-cell}
from astropy import time
from astropy import units as u
import numpy as np
from poliastro.bodies import Earth
from poliastro.core.propagation import func_twobody
from poliastro.examples import iss
from poliastro.plotting import OrbitPlotter3D
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import CowellPropagator
from poliastro.twobody.sampling import EpochsArray
from poliastro.util import norm
```
```{code-cell}
# More info: https://plotly.com/python/renderers/
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook_connected"
```
To provide an acceleration depending on an extra parameter, we can use **closures** like this one:
```{code-cell}
accel = 2e-5
```
```{code-cell}
def constant_accel_factory(accel):
def constant_accel(t0, u, k):
v = u[3:]
norm_v = (v[0] ** 2 + v[1] ** 2 + v[2] ** 2) ** 0.5
return accel * v / norm_v
return constant_accel
```
```{code-cell}
def f(t0, state, k):
du_kep = func_twobody(t0, state, k)
ax, ay, az = constant_accel_factory(accel)(t0, state, k)
du_ad = np.array([0, 0, 0, ax, ay, az])
return du_kep + du_ad
```
```{code-cell}
times = np.linspace(0, 10 * iss.period, 500)
times
```
```{code-cell}
ephem = iss.to_ephem(
EpochsArray(iss.epoch + times, method=CowellPropagator(rtol=1e-11, f=f)),
)
```
And we plot the results:
```{code-cell}
frame = OrbitPlotter3D()
frame.set_attractor(Earth)
frame.plot_ephem(ephem, label="ISS")
```
## Error checking
```{code-cell}
def state_to_vector(ss):
r, v = ss.rv()
x, y, z = r.to_value(u.km)
vx, vy, vz = v.to_value(u.km / u.s)
return np.array([x, y, z, vx, vy, vz])
```
```{code-cell}
k = Earth.k.to(u.km**3 / u.s**2).value
```
```{code-cell}
rtol = 1e-13
full_periods = 2
```
```{code-cell}
u0 = state_to_vector(iss)
tf = (2 * full_periods + 1) * iss.period / 2
u0, tf
```
```{code-cell}
iss_f_kep = iss.propagate(tf)
```
```{code-cell}
iss_f_num = iss.propagate(tf, method=CowellPropagator(rtol=rtol))
```
```{code-cell}
iss_f_num.r, iss_f_kep.r
```
```{code-cell}
assert np.allclose(iss_f_num.r, iss_f_kep.r, rtol=rtol, atol=1e-08 * u.km)
assert np.allclose(
iss_f_num.v, iss_f_kep.v, rtol=rtol, atol=1e-08 * u.km / u.s
)
```
```{code-cell}
assert np.allclose(iss_f_num.a, iss_f_kep.a, rtol=rtol, atol=1e-08 * u.km)
assert np.allclose(iss_f_num.ecc, iss_f_kep.ecc, rtol=rtol)
assert np.allclose(iss_f_num.inc, iss_f_kep.inc, rtol=rtol, atol=1e-08 * u.rad)
assert np.allclose(
iss_f_num.raan, iss_f_kep.raan, rtol=rtol, atol=1e-08 * u.rad
)
assert np.allclose(
iss_f_num.argp, iss_f_kep.argp, rtol=rtol, atol=1e-08 * u.rad
)
assert np.allclose(iss_f_num.nu, iss_f_kep.nu, rtol=rtol, atol=1e-08 * u.rad)
```
## Numerical validation
According to [Edelbaum, 1961], a coplanar, semimajor axis change with tangent thrust is defined by:
$$\frac{\operatorname{d}\!a}{a_0} = 2 \frac{F}{m V_0}\operatorname{d}\!t, \qquad \frac{\Delta{V}}{V_0} = \frac{1}{2} \frac{\Delta{a}}{a_0}$$
So let's create a new circular orbit and perform the necessary checks, assuming constant mass and thrust (i.e. constant acceleration):
```{code-cell}
orb = Orbit.circular(Earth, 500 << u.km)
tof = 20 * orb.period
ad = constant_accel_factory(1e-7)
def f(t0, state, k):
du_kep = func_twobody(t0, state, k)
ax, ay, az = ad(t0, state, k)
du_ad = np.array([0, 0, 0, ax, ay, az])
return du_kep + du_ad
orb_final = orb.propagate(tof, method=CowellPropagator(f=f))
```
```{code-cell}
da_a0 = (orb_final.a - orb.a) / orb.a
da_a0
```
```{code-cell}
dv_v0 = abs(norm(orb_final.v) - norm(orb.v)) / norm(orb.v)
2 * dv_v0
```
```{code-cell}
np.allclose(da_a0, 2 * dv_v0, rtol=1e-2)
```
This means **we successfully validated the model against an extremely simple orbit transfer with an approximate analytical solution**. Notice that the final eccentricity, as originally noticed by Edelbaum, is nonzero:
```{code-cell}
orb_final.ecc
```
## References
* [Edelbaum, 1961] "Propulsion requirements for controllable satellites"